function [] = haar_trainning(pos_dir,neg_dir,state_file,T)

% load state if exist
if exist(state_file, 'file') == 2
	state = load(state_file);
	haar_trainning_weight = state.haar_trainning_weight;
	haar_trainning_output = state.haar_trainning_output;
	POS = state.POS;
	NEG = state.NEG;
else
	[POS,NEG]  = load_trainning_data(pos_dir,neg_dir);

	% number of pos and neg trainning examples
	pos_num = size(POS,3)
	neg_num = size(NEG,3)

	haar_trainning_weight = [ones((pos_num),1).*(1/pos_num) ; ones((neg_num),1).*(1/neg_num)];
	haar_trainning_output = [];
	save(state_file);
end

for i = 1:T
	haar_trainning_weight = haar_trainning_weight ./ sum(haar_trainning_weight);

	[haar_trainning_weight current_best] = find_the_best_feature(haar_trainning_weight,POS,NEG);
	fprintf('iteration=%d,total=%d\n',i,T);

	haar_trainning_output(:,size(haar_trainning_output,2)+1)=current_best;

	save(state_file);
end

function [POS,NEG] = load_trainning_data(pos_dir,neg_dir)
pos_images = strcat(pos_dir,'/','*.pgm');
neg_images = strcat(neg_dir,'/','*.pgm');
pos_set = dir(pos_images);
neg_set = dir(neg_images);

num_pos = size(pos_set,1);
num_neg = size(neg_set,1);

for dir_i = 1:num_pos
        dir_name = pos_set(dir_i).name;
	path_name = strcat(pos_dir,'/',dir_name);
	img = double(imread(path_name));
	var_v = var(img(:));
	if var_v == 0
		var_v = 1;
	end
	img = (img - mean(img(:)))/sqrt(var_v);
	POS(:,:,dir_i)=cumsum(cumsum(img),2);
end

for dir_i = 1:num_neg
        dir_name = neg_set(dir_i).name;
	path_name = strcat(neg_dir,'/',dir_name);
	img = double(imread(path_name));
	var_v = var(img(:));
	if var_v == 0
		var_v = 1;
	end
	img = (img - mean(img(:)))/sqrt(var_v);
	NEG(:,:,dir_i)=cumsum(cumsum(img),2);
end



function [new_weight this_result]= find_the_best_feature(haar_trainning_weight,POS,NEG);

feature_type = [ 1,2; 2,1; 1,3; 3,1; 2,2];

HEIGHT = size(POS,1);
WIDTH = size(POS,2);

count = 0;

tranning_example_mask = [ones(size(POS,3),1);zeros(size(NEG,3),1)];
neg_total_weight = (~tranning_example_mask)'*haar_trainning_weight;
pos_total_weight = (tranning_example_mask)'*haar_trainning_weight;

err_best =1;
x_best = 0;
y_best = 0;
x_end_best=  0;
y_end_best=0;
feature_best=0;
detect_best=0;
p_best=0;
th_best=0;

for feature=1:size(feature_type,1)
	for x = 1:HEIGHT
		for y=1:WIDTH
			basic_x = feature_type(feature,1);
			basic_y = feature_type(feature,2);

			for explorer_x=x+basic_x-1:basic_x:HEIGHT
				for explorer_y=y+basic_y-1:basic_y:WIDTH
					[p,th,detect,error_v] = evalulate_this_feature(POS,NEG,x,y,explorer_x,explorer_y,feature, haar_trainning_weight, pos_total_weight, neg_total_weight);
					if error_v < err_best
						err_best = error_v;
						x_best = x;
						y_best = y;
						x_end_best = explorer_x;
						y_end_best = explorer_y;
						feature_best = feature;
						detect_best = detect;
						p_best = p;
						th_best=th;
						fprintf('new best detected: x=%d, y=%d, x_end=%d, y_end=%d, feature=%d(%d,%d), p=%d, th=%f -> err=%f\n',x_best,y_best,x_end_best,y_end_best,feature_best,basic_x,basic_y,p_best,th_best,err_best);
					end
					count = count + 1;
					if mod(count,10000) == 0
						fprintf('status update: count=%d, feature=%d\n',count,feature);
					end
				end
			end
		end
	end
end
count

wrong = xor(tranning_example_mask,detect_best);
right=~wrong;
beta = err_best/(1-err_best);
alpha = log(1/beta);
right = right * beta; 
right = right .* haar_trainning_weight;
wrong = wrong .* haar_trainning_weight;
haar_trainning_weight = right + wrong;
new_weight=haar_trainning_weight;
this_result=[x_best;y_best;x_end_best;y_end_best;feature_best;p_best;th_best;err_best;beta;alpha];

function [p,th,detect,error_v] = evalulate_this_feature(POS,NEG,x,y,explorer_x,explorer_y,feature, haar_trainning_weight, pos_total_weight,neg_total_weight)
num_pos = size(POS,3);
num_neg = size(NEG,3);

result_p = ones(num_pos,2);
result_n = zeros(num_neg,2);
result_p(:,2) = haar_evaluate_feature(POS,x,y,explorer_x,explorer_y,feature);
result_n(:,2) = haar_evaluate_feature(NEG,x,y,explorer_x,explorer_y,feature);

result = zeros(num_pos+num_neg,3);
result(:,1:2) = [result_p;result_n];
result(:,3) = haar_trainning_weight;
[~,indx]=sort(result(:,2));
result=result(indx,:);
%result=sortrows(result,2);

result_a = zeros(num_pos+num_neg,4);
result_a(:,1) = cumsum(result(:,1).*result(:,3));
result_a(:,2) = cumsum((~result(:,1)).*result(:,3));
result_a(:,3) = result_a(:,1)+neg_total_weight-result_a(:,2);
result_a(:,4) = result_a(:,2)+pos_total_weight-result_a(:,1);
[pos_min,pos_index] = min(result_a(:,4));
[neg_min,neg_index] = min(result_a(:,3));

detect = [result_p(:,2);result_n(:,2)];
if pos_min <= neg_min
	th = result(pos_index,2);
	p = 1;
	error_v = pos_min;
	detect = detect*p <= p* th;
else
	th = result(neg_index,2);
	p = -1;
	error_v = neg_min;
	detect = detect*p < p* th;
end
